R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter. # Tue Jun 7 19:35:18 2022 ——————————

source("tianfengRwrappers.R")

Warning in max(merge_expr$expr) : max里所有的参数都不存在；回覆-Inf

rat Frontier GSE174098 carotid

rat10x <- CreateSeuratObject(Read10X("./rat_scRNAseq/"), names.field = 2, names.delim = "-",
                                   project = "rat", min.cells = 10, min.features = 300) %>%
  PercentageFeatureSet(pattern = "^Mt-", col.name = "percent.mt") 

table(rat10x$orig.ident)
VlnPlot(rat10x,"nCount_RNA") /
VlnPlot(rat10x,"percent.mt") /
VlnPlot(rat10x, "nFeature_RNA")


rat10x <- rat10x %>% subset(subset = nFeature_RNA > 400 & nFeature_RNA < 4000 &
                              nCount_RNA > 1000 &  nCount_RNA < 30000 & percent.mt< 10) %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F) %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)

mouse10x <- CreateSeuratObject(Read10X("./celldiscovery_mouse_10x/1/"), names.field = 2, names.delim = "-",
                                   project = "mouse1", min.cells = 10, min.features = 300) %>%
  PercentageFeatureSet(pattern = "^mt-", col.name = "percent.mt") 

mouse10x2 <- CreateSeuratObject(Read10X("./celldiscovery_mouse_10x/2/"), names.field = 2, names.delim = "-",
                                   project = "mouse2", min.cells = 10, min.features = 300) %>%
  PercentageFeatureSet(pattern = "^mt-", col.name = "percent.mt") 

mouse10x <- merge(mouse10x,mouse10x2)

VlnPlot(mouse10x,"nCount_RNA") /
VlnPlot(mouse10x,"percent.mt") /
VlnPlot(mouse10x, "nFeature_RNA")

mouse10x <- mouse10x %>% subset(subset = nFeature_RNA > 500 & nFeature_RNA < 5000 &
                              nCount_RNA > 1000 &  nCount_RNA < 30000 & percent.mt < 10) %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F) %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)
# table(mouse10x$orig.ident)


umapplot(mouse10x)
f("Lmo2",mouse10x)

mouse coronary GSE131778

mouse_coronary_countmatrix <- read.csv("./GSE131776_mouse_scRNAseq.txt", sep = "\t")
func <- function(s) {
  paste0(strsplit(s, ".", fixed = T)[[1]][2], "_", strsplit(s, ".", fixed = T)[[1]][1])
}
colnames(mouse_coronary_countmatrix) <- lapply(colnames(mouse_coronary_countmatrix), func) # 拆分样本

mousecor <- CreateSeuratObject(counts = mouse_coronary_countmatrix,
                                   project = "mouse_cor", min.cells = 10, min.features = 300) %>%  PercentageFeatureSet(pattern = "^mt-", col.name = "percent.mt") 

# saveRDS(mousecor,"mousecor.rds")
table(mousecor$orig.ident)
VlnPlot(mousecor,"nCount_RNA") /
VlnPlot(mousecor,"percent.mt") /
VlnPlot(mousecor, "nFeature_RNA")


mousecor <- mousecor %>% subset(subset = nFeature_RNA > 400 & nFeature_RNA < 4000 &
                              nCount_RNA > 1000 &  nCount_RNA < 30000 & percent.mt < 10) %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F) %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)
table(mousecor$orig.ident)
saveRDS(mousecor,"mousecor.rds")

SMC2

mousecor <- readRDS("mousecor.rds")
f("Prdm16",mousecor, label.size = 7) + theme(legend.text = element_text(size = 20))

mouse carotid scRNAseq GSE155513

dataload

count_mats <- list.files("./GSE155513_RAW/")
count_mats <- count_mats[count_mats != "sampleinfo.txt"]
allList <- lapply(count_mats, function(file) {
  dd <- read.table(paste0("./GSE155513_RAW/", file), row.names = 1,stringsAsFactors = F)
  colnames(dd) <- as.character(dd['gene',])
  dd <- dd[-1,]
  CreateSeuratObject(
    counts = dd,
    project = file, min.cells = 10, min.features = 300
  )
})
# 合并seurat对象
mouse_carotid <- merge(allList[[1]], 
  y = allList[-1], add.cell.ids = count_mats,
  project = "mouse_carotid"
)
rm(allList)

# saveRDS(mouse_carotid,"mouse_carotid.rds")
# View(mouse_carotid)

process

SMC

mouse_carotid <- readRDS("mouse_carotid.rds")
umapplot(mouse_carotid)
umapplot(mouse_carotid,group.by = "orig.ident",label = F)
f("Prdm16",mouse_carotid) #SMC
f("Ly6a",mouse_carotid) #SEM-like cells

multi_featureplot(c("Bmp2","Bmp4","Bmp6"),mouse_carotid)

## BMP4 在这里EC的*大部分*中表达，而在人类样本中BMP4+ EC细胞是少数的
markers <- FindAllMarkers(mouse_carotid,logfc.threshold = 0.5,min.diff.pct = 0.2, only.pos = T)

stromal cells

mouse_carotid_stromal <- subset(mouse_carotid,idents = c(0,2,1,7))
mouse_carotid_stromal <- mouse_carotid_stromal %>% RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% FindClusters(resolution = 0.1)
mouse_carotid_stromal <- mouse_carotid_stromal %>% FindClusters(resolution = 0.2)

mouse_carotid_stromal <- readRDS("mouse_carotid_stromal.rds")
# saveRDS(mouse_carotid_stromal,"mouse_carotid_stromal.rds")

umapplot(mouse_carotid_stromal)
f("Dlx2",mouse_carotid_stromal) #Dlx2,Dlx5,Dlx6共同定位

mouse_carotid_stromal <- AddModuleScore(mouse_carotid_stromal, list(mmSMC2_marker), name = "SMC2_score")
f("SMC2_score1", mouse_carotid_stromal,min.cutoff = 0)

multi_featureplot(mmSMC2_marker[1:9],mouse_carotid_stromal,labels = NA, label = F)

分群表

table(group_tab[Dlx5poscells])/table(group_tab) 

GSM4705592_RPS003_matrix.txt.gz GSM4705593_RPS004_matrix.txt.gz GSM4705594_RPS011_matrix.txt.gz GSM4705595_RPS012_matrix.txt.gz GSM4705596_RPS007_matrix.txt.gz 
                    0.080645161                     0.076949502                     0.130000000                     0.029882604                     0.058282209 
GSM4705597_RPS008_matrix.txt.gz GSM4705598_RPS001_matrix.txt.gz GSM4705599_RPS002_matrix.txt.gz GSM4705600_RPS017_matrix.txt.gz GSM4705601_RPS018_matrix.txt.gz 
                    0.004746835                     0.039230575                     0.003976143                     0.093492209                     0.095634096 
GSM4705602_RPS013_matrix.txt.gz GSM4705603_RPS014_matrix.txt.gz GSM4705604_RPS015_matrix.txt.gz GSM4705605_RPS016_matrix.txt.gz 
                    0.066773504                     0.019255456                     0.042503503                     0.030303030 

图

ggsave("./fig7_mouse/cir_mouse_carotid_stromal2.png" ,plot = p, height = 12, width = 16,device = png)

Warning: Removed 1 rows containing missing values (geom_text).

human ds2

p <- multi_featureplot(c("FRZB","SOST","DLX5","DLX6"), ds2, labels = NA, label.size = 6)
ggsave("refds2_SMC2_carotid_stromal.png",plot = p, height = 7, width = 7,device = png)

Warning: Removed 1 rows containing missing values (geom_text).

human bulk RNA-seq GSE120521 carotid stable/unstable FPKM

fpkm2tpm <- function(fpkm){
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}

fpkm_matrix <- read.csv("GSE120521_FPKM.csv")
# fpkm_matrix <- distinct(fpkm_matrix) #去除重复行
fpkm_matrix <- fpkm_matrix[!duplicated(fpkm_matrix$name),]
rownames(fpkm_matrix) <- fpkm_matrix$name
fpkm_matrix$name <- NULL

tpm_matrix <- apply(fpkm_matrix, 2, fpkm2tpm)
colSums(tpm_matrix)

group_file <- c("stable","unstable","stable","unstable",
                "stable","unstable","stable","unstable")
boxplot(tpm_matrix, las = 2)

expr_mat <- tpm_matrix[!apply(tpm_matrix, 1, function(x){sum(floor(x) == 0)>3}),]

boxplot(expr_mat, las = 2)

library(limma)
expr_mat <- normalizeBetweenArrays(expr_mat)
expr_mat <- log2(expr_mat+1) #使用log2 scale

#PCA
library(ggfortify) 
df <- as.data.frame(t(expr_mat)) 
df$group <- group_file 
autoplot(prcomp(df[,1:(ncol(df)-1)]), data=df, colour = 'group')+ theme_bw() 

# If the sequencing depth is reasonably consistent across the RNA samples, then the simplest and most robust approach to differential exis to use limma-trend.

fit <- lmFit(expr_mat, group_file)
fit <- treat(fit, lfc=log2(1.2), trend=TRUE)
topTreat(fit, coef=ncol(design))

library(ggpubr)
dat <- expr_mat
design <- model.matrix(~factor(group_file))
fit <- lmFit(dat, design)
fit <- eBayes(fit)
# options(digits = 4)
topTable(fit,coef=2,adjust='BH')
deg <- topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg) 

write.csv(deg,"./datatable/stable vs unstable.csv")
## look up FRZB, SOST, PRDM6, OGN
ds1markers[ds1markers$cluster == "SMC2",]$gene

ds2markers[ds2markers$cluster == "SMC3",]$gene

smc2markers <- intersect(ds1markers[ds1markers$cluster == "SMC2",]$gene, ds2markers[ds2markers$cluster == "SMC3",]$gene)


deg[intersect(smc2markers, rownames(deg)),]

ds1markers[ds1markers$cluster == "SMC1",]$gene

ds2markers[ds2markers$cluster == "SMC1",]$gene

SMC1markers <- intersect(ds1markers[ds1markers$cluster == "SMC1",]$gene, ds2markers[ds2markers$cluster == "SMC1",]$gene)

deg[intersect(SMC1markers, rownames(deg)),]

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.